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We theoretically demonstrate the existence of frozen light in near-zero-index media with cubic 
nonlinearity. Light is stopped to a standstill owing to the divergent wavelength and the vanish¬ 
ing group velocity, effectively rendering, through nonlinearity, a positive-epsilon trapping cavity 
carved in an otherwise slightly-negative-epsilon medium. By numerically solving Maxwell’s equa¬ 
tions, we find a soliton-like family of still azimuthal doughnuts, which we further study through an 
adiabatic perturbative theory that describes soliton evaporation in lossy media or condensation in 
actively pumped materials. Our results suggest applications in optical data processing and storage, 
quantum optical memories, and soliton-based lasers without cavities. Additionally, near-zero-index 
conditions can also be found in the interplanetary medium and in the atmosphere, where we provide 
an alternative explanation to the rare phenomenon of ball-lightning. 

PACS numbers: 42.65.Tg, 42.65.Wi, 42.79.Gn, 71.45.Lr 


Introduction - Optical beams are generally unbound 
in free-space and in bulk media. Guidance and full spa¬ 
tial confinement of light are usually achieved by means 
of waveguides, mirrors, resonators, and photonic crys¬ 
tals. Alternatively, nonlinear self-organization can be ex¬ 
ploited to compensate for diffraction of optical beams or 
dispersive broadening of pulses, enabling the formation of 
spatial and temporal solitons, respectively W- Spatial 
self-trapping occurs in several optical systems, including 
photorefractive media 13 E], liquid crystals 0 m, and 
metamaterials um- Remarkably, nonlinearity can act 
simultaneously on temporal and spatial domains to com¬ 
pensate for both diffraction and dispersion, thus enabling 
the formation of light bullets, spatio-temporal dough¬ 
nuts, and X-shaped waves P3HI3- 

Physical systems enabling either slow or fast light 
H8| naturally enhance radiation-matter interaction, thus 
boosting nonlinear processes that can be efficiently used 
for active light control m , all-optical switching, and 
modulation [2Qll2Tj . In particular, near-zero-index (NZI) 
media can slow down light propagation [22, [23] and en¬ 
able extreme nonlinear dynamics [25] . enhanced second 
and third harmonic generation ESI, active control of tun¬ 
neling optical switching, and bistable response [28] . 
These materials naturally exist in nature, for example 
plasmas, transparent conductors, and metals near the 
their bulk plasma frequency uo p [29] . Besides, they can 
be artificially realized as waveguides close to modal cut¬ 
off [30] . using surface phonon polaritons in GaAs quan¬ 
tum wells m , or by engineering sub wavelength metal¬ 
lic nanowires, nano-spheres, or nano-circuits embedded 
in dielectric matrices. The latter strategy has enabled 
the development of epsilon-near-zero (ENZ) metamateri¬ 
als, which have been investigated for applications such as 
enhanced transmission [32] , cloaking [33] , energy squeez¬ 
ing in narrow channels [34] . and subwavelength imaging 
[35] 3b. The ENZ regime is inevitably associated with 
high dispersion and is therefore accompanied by absorp¬ 


tion, which can be suppressed by embedding externally 
pumped active inclusions in the NZI medium 123. 

In this Letter, we investigate self-organization of light 
in NZI media with Kerr-like instantaneous nonlinearity. 
In particular, we reveal the existence of fully confined 
doughnut-shaped solitons with vanishing Poynting vector 
and angular momentum. In practice nonlinearity enables 
digging a three-dimensional cavity for light, which in turn 
remains frozen and self-trapped. We study the effect of 
loss on stationary light doughnuts by developing a fully 
numerical soliton perturbative theory, finding that they 
evaporate over time due to inelastic absorption: their am¬ 
plitude decreases, their frequency blueshifts slightly, and 
their radius increases. Conversely, if externally pumped 
active inclusions with inversion of population are embed¬ 
ded within the NZI medium, the opposite scenario takes 
place and azimuthal doughnuts condensate over time. 
These findings demonstrate the possibility to freeze light 
beams in ENZ media, with potential applications in op¬ 
tical data processing and storage, quantum optical mem¬ 
ories, and NZI lasers operating without cavities. While 
confinement in random lasers is generally brought by An¬ 
derson localization [38H4Q] , in the case of NZI media self¬ 
trapping is provided by optical nonlinearity in the form of 
solitary waves. Interestingly, ENZ conditions are found 
also in the interplanetary medium and in the atmosphere, 
and we argue that our theoretical results may provide in¬ 
sight into ball-lightning (BL) formation [4TH43] . 

Model - We consider a generic NZI medium with 
Drude temporal response and instantaneous Kerr-like 
nonlinearity. Both of these ingredients ensue from free- 
particle temporal dynamics, which is characteristic of 
plasmas, metals, transparent conductors, and ENZ meta¬ 
materials, all examples of NZI media. In particular, Kerr- 
like nonlinearity naturally arises from the ponderomotive 
force in plasmas and metals [44] . and is well represented 
by the constitutive relation between the displacement 
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FIG. 1: (Color online) Dispersion relation k(uS) (cyan right 
y- axis) and phase and group velocities (v/ and v g , red left 
y- axis) of TEM waves in the linear loss-less limit (x 3,7 — > 0). 
All quantities are plotted in dimensionless units: the angular 
frequency uj is normalized to the plasma frequency uj p , the 
wave-vector k is normalized to k p = uj p /c , and the phase 
and group velocities are normalized to the speed of light in 
vacuum c. The red dashed line indicates the dispersion-less 
limit v f = v g = c. 


vector R e[D(t)\ and the electric field Re [£(£)]: 

f'OG 

V{t) = eo / e(t')£(t - t')dt' + (1) 

Jo 

where eo is the vacuum permittivity, X 3 is the nonlinear 
susceptibility of the medium, e(r) = S(r)-\-uJ p (l—e~ ryT )/ / y 
is the Drude temporal response function, S(r) is the Dirac 
delta-function, uj p is the plasma frequency, and 7 is the 
temporal damping rate due to inelastic collisions. Optical 
propagation is governed by the wave equation 

VxVxf = -/i 0 <9 t 2 X>, (2) 

where fi 0 is the vacuum permeability. In the linear limit 
X 3 —x 0, homogeneous transverse electromagnetic (TEM) 
waves are solutions of Eq. Q given by £ = eoe 2k r_zu;t , 
where k • eo = 0. The angular frequency uj and the 
wave-vector k satisfy the dispersion relation k(uj) = 
(uj/c )Ei uj), where c is the speed of light in free space and 
e(uj) = 1 — ojp/uj(oj + 27 ) is the frequency-dependent di¬ 
electric constant, which is given by the Fourier transform 
of the Drude temporal response function e(r). The linear 
dispersion relation of TEM waves k{uS) is depicted in Fig. 
[I] in the lossless limit 7 —x 0 , together with the phase and 
group velocities Vf(uj) = uj/k(uj), v g (uj) = duj/dk. Note 
the cutoff of TEM waves at the plasma frequency uj = u; p , 
where the medium enters the ENZ regime, the phase ve¬ 
locity diverges, and the group velocity vanishes [221123] , 
Homogeneous nonlinear modes - Owing to the van¬ 
ishing group velocity, nonlinear effects are dramatically 



FIG. 2: (Color online) (a) Nonlinear dispersion of zero-index 
homogeneous modes existing for angular frequencies uj smaller 
than the cutoff uj p . The electric field amplitude E 0 is normal¬ 
ized to the scaling electric field Es = . (b,c) Maximum 

instability growth a' M (normalized to the plasma frequency 
cj p ) as a function of the perturbing wave-vector modulus q 
(normalized to k p — uj p /c) for several directions in the recip¬ 
rocal space: (b) 0 — 7r/4, (f> = 0,7r/4, 7t/2, and (c) <p = 0, 
and 6 — 0,7r/4, 7t/2. (d) Contour-plot of the maximum insta¬ 
bility growth a M (normalized to the plasma frequency cj p ) as 
a function of 0, 4> for a fixed perturbing wave-vector modulus 
q/kp = 0.1. 


enhanced in the ENZ regime [25, 26.. For uj < cj p , ho¬ 
mogeneous modes with vanishing wavenumber, infinite 
phase velocity, and zero group velocity can be found 
by neglecting damping and setting £ = E 0 e~ luJt ^ with 
Eo = 2c(uj) /(3x3)• The resulting dispersion rela¬ 
tion is plotted in Fig. |2|a). We find zero-index ho¬ 
mogeneous modes to have a cutoff at the plasma fre¬ 
quency cj p , where the electric field amplitude drops to 
zero. In order to evaluate the stability of homoge¬ 
neous modes, we perturb them with small-amplitude 
waves: £ = [E 0 + JEie iq,r+at + e~ iujt , 

where JEi, 5 E 2 are the perturbation amplitudes with 
wave-vector q and temporal growth eigenvalue a. In¬ 
serting this expression in Eqs. 0 and <0 and retaining 
only the lowest-order terms in S Ei and c)E 2 , we find a ho¬ 
mogeneous system of linear equations, whose non-trivial 
solutions are signaled by the vanishing of the secular de¬ 
terminant [24 . This condition determines the complex 
temporal eigenvalues a. Instabilities are then associated 
with positive real parts of the eigenvalue a , indicating un¬ 
bound amplification of the perturbation. We plot results 
of the stability analysis in Figs |2jb-d) , and in particular, 
we depict the maximum of the real part of the eigen¬ 
value, I n analogy to standard modulation instabil¬ 
ity in ID paraxial systems [3], the gain spectrum of the 
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perturbations is non-vanishing within a finite wavevec- 
tor window and is peaked at a characteristic wavevector 
modulus. However, in contrast to ID paraxial systems, 
the gain spectrum is 3D and has a non-trivial dependence 
on polar and azimuthal angles (0,0) of the perturbation 
wavevector q. 

Still azimuthal doughnuts - The modulation instabil¬ 
ity scenario strongly suggests the presence of still 3D soli- 
tons in NZI media. In order to verify this hypothesis, we 
transform Eq. § into spherical coordinates and search 
for azimuthally-polarized solutions: £ = F^(r, 0)e _zu;t 0- 
As Eq. © is invariant under a constant phase shift, 
without any loss of generality we can assume that the 
electric field envelope is real F^(r, 0) G meaning that 
we are seeking non-propagating solutions which are not 
accompanied by a phase flow. Indeed, assuming that 
such solutions exist, we show that the Poynting vector 
vanishes thoroughly (see SI EH). Besides, we seek lo¬ 
calized soliton-like solutions vanishing at r —>■ oo and 
atr = 0 , 0 = 0 , 7 r owing to the azimuthal polarization. 
Upon examination of the asymptotical expansion of Eq. 
© for r —>> ex), we find that 3D soliton-like azimuthal 
solutions can actually exist only in the ENZ regime (see 
SI [24]). Thus, we discretize derivatives with respect to 
the radius r n and the polar angle 0 m and then transform 
the differential wave equation for the electric field into 
a nonlinear algebraic system for the electric field ampli¬ 
tudes E^ nrn in the two-dimensional grid r n ,0 m (see SI 
EH). We solve this nonlinear algebraic system by means 
of an iterative Newton-Raphson algorithm, and find a 
family of still azimuthal doughnuts [see Fig. ©a)] for 
uj < cj p , which presents a cutoff at cj p , where the soli- 
ton loses localization and its amplitude vanishes. The 
frequency-dependent maximum amplitude and the corre¬ 
sponding radius of the still doughnut family are plotted 
in Fig. [3jb), while a r-6 contour-plot of the squared elec¬ 
tric field profile |F 0 (r)/Fs | 2 (normalized to the scaling 
field Fg = 1/y/xs) of the still doughnut at uj/uj p = 0.995 
is depicted in Fig. ^c). The total dielectric permittiv¬ 
ity profile eT(cJ,r) = e(oj) + (3/2)x3\E^(cj, r )| 2 is shown 
in Fig. [3|d). Importantly, in the soliton existence do¬ 
main u < cu p , the linear dielectric constant is negative 
e(ca) < 0 , and thus, at long radius where the electric 
field amplitude is small, the NZI medium is metal-like. 
Conversely, in the volume around the radius r max for 
which the electric field is maximum, nonlinearity is non- 
negligible and the total dielectric permittivity is positive 
ex(r m ax) > 0 (dielectric-like). From here we see that the 
existence of still azimuthal doughnuts originates in the 
extraordinary ability of nonlinearity to dig a dielectric¬ 
like 3D cavity within a metal-like environment. This sce¬ 
nario is unique of NZI media, which prevent propagation 
of the fields outside the induced-dielectric trapping cav¬ 
ity. We emphasize that modulation instability enables 
the excitation of non-propagating solitons starting di¬ 
rectly from unstable homogeneous waves with frequency 



FIG. 3: (Color online) (a) Iso-surface |F^(r)/Fs| 2 = 

0.005 of a still doughnut with maximum squared amplitude 
|Fm/Fs | 2 = 0.01, where Es = is the scaling field am¬ 

plitude, excited at an angular frequency uj/uj p — 0.9992. (b) 
Soliton maximum amplitude Fm/Fs (red left y- axis) and its 
corresponding radius r ma x (cyan right y- axis) as a function of 
angular frequency cj/ujp. (c) Contour-plot in the r-6 plane 
of the dimensionless intensity profile |F^(r, 0)/Fs| 2 and (d) 
total dielectric permittivity profile ct (r*, 6) (red surface) asso¬ 
ciated with the still doughnut of (a). The blue plane in (d) 
represents the metal-dielectric transition plane ct (r*, 0) = 0. 
All quantities are plotted in dimensionless units: the angu¬ 
lar frequency uj is normalized to the plasma frequency cj p , 
while spatial coordinates are normalized to the inverse of the 
plasma wave-vector k ~ 1 . 


falling in the ENZ regime. 

Doughnut evaporization/condensation - In standard 
transparent media, the main quantity accounting for op¬ 
tical propagation is the Poynting vector, representing 
the photon flux per unit time. For our trapped soli¬ 
tons, the Poynting vector is thoroughly vanishing (see 
SI [24]), so we describe doughnut self-trapping through 
the optical-cycle-averaged density of electromagnetic en¬ 
ergy u = (l/ 2 )£ • V. Now, if absorption is taken into 
account, the energy density is expected to be damped 
and vanish exponentially over time. A numerical ver¬ 
ification of this hypothesis could consist in temporally 
evolving Eq. © with the doughnut initial condition. 
However, temporal evolution requires nonlinear 3D finite- 
difference-time-domain (FDTD) numerical simulations, 
which are computationally demanding. Besides, tradi¬ 
tional approaches used in dielectric and plasmonic waveg¬ 
uides mm relying on the slowly-varying-envelope ap¬ 
proximation (SVEA) can not be used, as the SVEA does 
not hold in the ENZ regime [25]. Instead, we have devel¬ 
oped a soliton perturbation theory (see SI [24]) capable of 
accounting for both damping and amplification (e.g., in 
systems containing externally pumped active inclusions 
within the NZI medium) under the assumption that (i) 
damping (7 > 0 ) or (ii) gain (7 < 0 ) are much smaller 
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FIG. 4: (Color online) (a) Soliton maximum amplitude 
Em/E s (red left y- axis) and its corresponding dimensionless 
radius k p r m ax (cyan right y- axis) as a function of dimension¬ 
less time where 7 is a phenomenological damping. (b-d) 
Iso-surfaces \E c f ) (r)/Es\ 2 = 0.001 of the time-evolving dough¬ 
nut with initial condition as in Fig. [3ja) for (b) 7 1 = 0 , (c) 
jt = 1, and (d) 7 t — 2 . 


than the soliton angular frequency uj. We further as¬ 
sume that the temporal evolution of the still doughnut 
adiabatically follows the soliton family, finding that the 
soliton amplitude (i) decays or (ii) increases over time fol¬ 
lowing the exponential law Eyiif) = Foe -7 */ 2 , where Eq 
is the initial field amplitude and 7 is a phenomenologi¬ 
cal absorption/pumping rate. Accordingly, the doughnut 

(i) expands and blueshifts or (ii) shrinks and redshifts 
in either case (see SI [24] ). The time-dependent field 
amplitude (blue left y- axis) and doughnut radius (red 
right y- axis) are plotted in Fig. |4^a) for a representative 
example, along with three snap-shots of the iso-surface 
|F^(r)/F?s| 2 = 1 x 10 -3 at different times in Figs. EF-d), 
where we have assumed as initial condition the doughnut 
of Fig. |3ja) and (i) damping (7 > 0) (For the full tem¬ 
poral evolution see movie in the SI [24] ). The doughnut 
evaporates over time, as its amplitude decreases and its 
radius increases. The gain scenario (ii) (7 < 0) can be 
interpreted by inverting the temporal direction, so that 
the doughnut condensates over time, as its amplitude in¬ 
creases and its radius decreases. 

Ball-lightnings? - BLs are rare lightning events with 
hitherto unknown theoretical explanation [4TH43] . BLs 
emit broadband radiation and can either propagate or 
stand still. Initially considered as myth, BLs have puz¬ 
zled scientists for centuries and their existence has been 
questioned until the first recent experiment able to mea¬ 
sure their spectrum [43]. Understanding of the nature of 
BLs is still unsatisfactory as they can not be easily repro¬ 
duced in laboratory. Among the several theories trying 
to explain their nature, the so-called maser-caviton the¬ 
ory m suggests that BLs are localized high-field soli- 
tons forming a cavity surrounded by plasma. Indeed, 


during thunderstorms, atmosphere can get ionized and 
become a NZI medium with a plasma frequency falling 
in the terahertz-microwave spectral region, where rota¬ 
tional levels of water can be excited. The ensuing emit¬ 
ted radiation is thought to remain self-trapped and heat 
up the air, thus emitting broadband blackbody radiation 
m- This theory explains some general aspects of BLs, 
but it does not provide any quantitative description of 
the self-induced soliton cavity. Following our rigorous 
calculations, we speculate that BLs may ensue from a 
self-organization process in the ENZ regime, where we 
observe the existence of still doughnut solitons, as dis¬ 
cussed above. The actual spherical shape of BLs observed 
in experiments [43] may be due to mixed polarization, 
higher order nonlinear effects, or the intrinsically inco¬ 
herent nature of radiation emitted in the atmosphere. 
The ENZ condition would explain the infrequency of the 
phenomenon and provides an insightful signature for ex¬ 
perimental investigations. 

Conclusions - Our investigation of self-organization 
phenomena in NZI media with cubic nonlinearity has 
resulted in the demonstration that zero-index nonlinear 
waves are unstable in all spatial directions and that still 
azimuthally polarized self-trapped doughnuts can be ex¬ 
cited. We have discussed the existence domain of this 3D 
soliton family with thoroughly vanishing Poynting vector 
and provided details on its characteristics. Besides, we 
have studied the effect of loss/amplification, finding that 
still light doughnuts evaporate/condensate over time, re¬ 
spectively. Our model applies to any NZI medium with 
cubic nonlinearity and our results are universal as they 
are rescaled to the relevant physical quantities (plasma 
frequency cj p , plasma wave-vector fc p , Kerr coefficient X 3 ) 
of any specific medium in this regime (e.g., metals, trans¬ 
parent conductors, plasmas, and metal-dielectric ENZ 
metamaterials). Our findings pave the way for the devel¬ 
opment of novel applications in optical data processing 
and storage, the realization of quantum optical memo¬ 
ries, and the design of soliton-based lasers without cav¬ 
ities. Incidentally, NZI conditions can be found also in 
the interplanetary medium and in the atmosphere, and 
we have discussed possible relationships between our re¬ 
sults and ball-lightning formation. 
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